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ABSTRACT 
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dynamics  of  the  observer  may  have  a  significant  effect  on  the  computed  gain  margin 
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I.  INTRODUCTION 

Accurate  path  control  of  surface  ships  and  underwater  vehicles  along  prescribed 
geographical  paths  is  a  basic  problem  that  is  becoming  increasingly  important,  par- 
ticularly as  the  missions  of  ocean  vehicles  become  more  complicated  with  strict 
requirements  for  performance.  In  order  for  a  control  law  to  be  able  to  perform  its 
mission  in  a  realistic  operational  scenario  it  has  to  be  robust  enough  so  that  it  can 
maintain  stability  and  accuracy  of  operations  in  the  presence  of  modeling  errors  and 
environmental  uncertainties.  The  robustness  properties  of  the  design  are  particu- 
larly important  due  to  the  unpredictable  nature  of  the  ocean  environment  and  the 
changes  in  the  hydrodynamic  characteristics  of  the  vehicle  during  turning,  changes 
in  the  forward  speed,  or  operations  in  proximity  to  other  objects  in  the  area.  For 
these  reasons,  there  exists  a  need  for  the  analysis  of  the  robustness  characteristics 
of  a  particular  control  law  design  and  the  establishment  of  a  rational  operational 
envelope  based  on  stability  and  performance  criteria.  Previous  studies  [Parsons 
and  Cuong  (1977)]  showed  that  gain  adaption  is  highly  desirable  due  to  changes 
in  the  linearized  vehicle  hydrodynamics  with  different  operation  conditions,  such 
as  depth  under  keel.  The  resulting  adaptation  scheme  [Parsons  and  Cuong  (1980)] 
required  significant  vehicle  motion,  which  may  be  undesirable  when  operating  in 
restricted  waters,  or  in  object  recognition  and  localization  tasks.  Integral  control 
techniques  [Parsons  and  Cuong  (1981)]  proved  quite  effective,  but  neglected  the 
nonlinear  behavior  of  the  vehicle,  which  becomes  very  important  at  low  speeds  and 
hover  operations.  Model  based  compensators  exhibit  robust  behavior  under  con- 
ditions of  parameter  uncertainty,  which  is  as  good  as  the  classical  linear  quadratic 


regulators  for  linear  output  feedback  systems  [Healey  (1992)].  Alternatively,  sliding 
mode  controllers  exhibit  very  robust  characteristics  given  an  estimate  of  the  param- 
eter uncertainty  and/or  disturbances  [Papoulias  and  Healey  (1992)],  [Yoerger  and 
Slotine  (1985)].  Sliding  mode  control,  however,  does  not  offer  an  infinitely  robust 
design  and  it  suffers  from  a  series  of  bifurcation  phenomena  and  loss  of  stability 
unless  proper  care  is  exercised  [Papoulias  (1991)]. 

In  this  work  we  analyze  the  problem  of  loss  of  stability  of  a  marine  vehicle 
under  cross  track  error  control  in  the  presence  of  mathematical  versus  actual  system 
mismatch.  For  demonstration  purposes,  variations  in  the  heading  angle  control 
gain  are  studied.  Previous  studies  [Oral  (1993)]  concentrated  on  system  response 
assuming  perfect  and  complete  state  measurement.  Particular  emphasis  in  this  work 
is  placed  on  analyzing  the  effects  of  observer  design  on  system  response  after  initial 
loss  of  stability  of  straight  line  motion.  The  main  loss  of  stability  cases  analyzed 
here  occur  in  the  form  of  generic  bifurcations  to  periodic  solutions  [Guckenheimer 
and  Holmes  (1983)].  We  use  center  manifold  reduction  techniques  and  averaging 
in  order  to  capture  the  stability  properties  of  the  resulting  limit  cycles  [Chow  and 
Mallet-Paret  (1977)].  It  is  shown  that  the  dynamics  of  the  observer  may  have  a 
significant  effect  on  the  computed  gain  margin  of  the  control  system  depending 
on  the  particular  basis  used.  All  computations  in  this  work  are  conducted  for  the 
NPS  autonomous  underwater  vehicle  [Bahrke  (1992)]  and  all  results  are  presented  in 
standard  dimensionless  quantities  with  respect  to  vehicle  length,  7.3  ft,  and  nominal 
forward  speed,  2  ft/sec. 


II.  PROBLEM  FORMULATION 

A.     EQUATIONS  OF  MOTION 

The  linear  maneuvering  equations  of  motion  of  a  marine  vehicle  in  the  hori- 
zontal plane  are  written  in  dimensionless  form  as, 

m(v  +  r  +  xGr)  =  Y^r  +  Y^v  +  Yrr  +  Yvv  +  Y68  (2.1) 

I2r  +  mxG[y  +  r)  =  N*r  +  Nbv  +  Nrr  +  Nvv  +  NSS  (2.2) 

where  all  symbols  are  explained  in  the  nomenclature.  Equations  (2.1)  and  (2.2)  can 
be  used  to  derive  a  second  order  transfer  function  between  the  rudder  angle  6  and 
yaw  rate  r.  For  low  frequency  maneuvering  motions  this  second  order  equation  can 
be  approximated  by  expanding  in  Taylor  series  and  keeping  the  first  order  terms 
only.  The  result  is 

r  =  ar  +  bS  (2.3) 

Equation  (2.3),  which  is  sometimes  referred  to  as  Nomoto's  first  order  model,  is 
particularly  useful  in  control  system  design  since  no  sway  velocity  feedback  is  neces- 
sary. This  equation  predicts  linear  variation  of  the  steady  state  turning  rate  versus 
rudder  angle.  In  reality,  the  r-6  curve  has  characteristics  of  softening  spring  mainly 
due  to  speed  loss  during  turning.  To  account  for  this  a  modified  version  of  equation 
(2.3)  is  used, 

r  =  ar  +  a3r3  +  b6  (2.4) 

where  a3  is  usually  determined  from  steady  state  results.  Finally,  the  model  is 
complete  by  the  incorporation  of  the  kinematic  equations, 

^  =  r  (2.5) 
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y  =  sin  tf  (2.6) 

where  \P  is  the  vehicle  heading,  and  y  is  the  cross  track  error  off  a  desired  straight 
line  path. 

B.     COMPENSATOR  DESIGN 

In  control  theory  it  is  known  that  the  eigenvalues  of  the  controller  are  not 
affected  by  the  eigenvalues  of  the  observer.  This  allows  us  to  design  the  controller 
and  observer  separately  which  is  known  as  the  separation  principle.  The  combination 
is  called  a  compensator. 

Equations  (2.3),  (2.5),  and  (2.6)  govern  the  steering  control  of  the  model  used 
in  this  section.  The  control  law  can  be  expressed  as, 

6  =  6MLtUnh(-^-)  (2.7) 

where  around  the  nominal  state  vp  =  r  =  y  =  0we  have 

60  =  K*V  +  Krr  +  Kyy  (2.8) 

6  is  the  rudder  angle  and  A*,  Kr,  and  Ky  are  the  control  gains  of  the  system. 
The  linear  control  law  is  60.  The  rudder  angle  6  is  defined  by  a  hyperbolic  tangent 
function  to  include  the  saturation  to  our  problem  as  shown  in  Figure  2.1.  Saturation 
occurs  at  ^sat,  which  is  the  saturation  limit  generally  taken  as  0.4  rad. 

The  linearized  form  of  equations  of  motions  in  the  vicinity  of^  =  r  =  y  =  0 
are, 

$  =  r  (2.9) 

r  =  ar  +  b60  (2.10) 

y  =  V  (2.11) 


A   5 

5 

ca: 

y 

«. 

-  —  o 

sat 

Figure  2.1:     Saturation  in  6. 


These  equations  can  be  expressed  in  state  space  form  as 


X  =  AX  +  Bu 


where 


X  = 


r 

y  . 


is  the  state  vector, 


A  = 


0 

1    0  1 

0 

a    0 

1 

0    0 

is  the  open  loop  dynamics  matrix  and 

B  = 

is  the  control  distribution  vector. 
The  observer  equations  are 


X  =  AX  +  Bu  +  L(Y  -  CX) 


(2.12) 


(2.13) 


where  X  is  the  estimate  of  X,  Y  is  the  output  of  the  system  Y  =  y,  and  C  is  the 
sensor  vector  C  =  [0  0  1]. 


The  error  in  the  estimate  of  X  is  defined  by 


X  =  X  -X 


Using  equations  (2.12),  (2.13),  and  (2.14)  we  can  obtain 


X  =  (A  -  LC)X 


We  can  rewrite  equations  (2.9),  (2.10),  and  (2.11)  in  the  form  of 
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The  observer  gains  are, 


L  = 


Cif/ 
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After  performing  the  matrix  operations  we  obtain 
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Using  equation  (2.13)  we  can  rewrite  equation  (2.8)  as  follows, 


60  =  K*(*  -  *)  +  Kr(r  -  r)  +  Ky{y  -  y) 


Finally,  we  can  write  our  compensator  equations  in  the  form 
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(2.15) 


(2.16) 
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If  we  look  at  the  matrix  carefully  we  will  see  that  it  is  in  the  form 

'  X'       \  A-BK       B 

k\~[      o        A- 

which  has  the  following  characteristic  equation, 

det[A  -  BK  -  si]  det[A  -  LC  -  al]  =  0 

This  indicates  that  the  dynamics  of  the  observer  are  completely  independent  of  the 
dynamics  (eigenvalues)  of  the  controller.  Thus  K  and  L  can  be  designed  separately. 

C.     CALCULATION  OF  CONTROL  GAINS 


A  is  the  Jacobian  matrix  of  the  system 

A  = 
The  characteristic  equation  of  the  matrix  A  is 


0  1  0 
bKv    a  +  bKr    bKy 

1  0  0 


A3  -  (a  +  bKr)\2  -  blUX  -  bKy  =  0 
If  the  desired  characteristic  equation  has  the  general  form 


A3  +  a2A2  +  ax\  +  a0  =  0 


the  control  gains  can  be  found  as 


K9    -    -T 


Kr    = 


a2  +  a 


Ky    -    -- 

The  desired  characteristic  equation  can  be  written  with  respect  to  the  desired  natural 
frequency  and  some  optimum  coefficients.  The  ITAE  criterion  for  a  third  order 
equation  is 

s3  +  1.75iyns2  +  2A5wls  +  wzn  =  0 


where  wn  is  the  desired  controller  natural  frequency. 

Therefore  the  control  gains  can  be  calculated  for  a  given  natural  frequency,  as 

ax    =    2.l5wl 
c*2    =    1.75u;n 

ao    =    wl 

D.     CALCULATION  OF  OBSERVER  GAINS 

If  we  define  A  as  the  Jacobian  matrix  of  the  system 

0    1     -£9 


A  = 


0  a    -£r 

1  0     -I 


v   J 


the  characteristic  equation  of  the  matrix  A  is 

A3  +  {£y  -  a)  +  (£*  -  a£y)X  +  (£r  -  a£*)  =  0 
If  the  desired  characteristic  equation  has  the  general  form 

A3  +  72A2  +  7iA  +  7o  =  0 
the  observer  gains  can  be  found  as 

£y     =     a +  72 
£y    =    a£y  +  7! 

£r    =    a£<i,  +  70 

Applying  the  ITAE  criteria,  observer  gains  can  be  calculated  for  a  given  natural 
frequency  as 

7i    =    2.15^0 

72    =    1.75u;no 

7o    =    w3n0 
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where  wno  is  the  observer  natural  frequency. 

E.     CHARACTERISTICS  OF  ITAE  CRITERIA 

In  the  calculations  of  gains  we  applied  the  ITAE  criteria.  If  we  look  at  Figure 
2.2  for  the  step  response  of  ITAE,  we  see  that  the  response  gets  faster  as  the  natural 
frequency  increases.  For  example,  the  settling  time  is  10  normalized  seconds,  or  10 
seconds  for  wn  =  1,  1  second  for  wn  =  10,  and  so  on. 


STEP  RESPONSE  of  ITAE 
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Figure  2.2:     Step  response  of  ITAE. 
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III.  HOPF  BIFURCATION 

A.     INTRODUCTION 

An  important  quantity  in  assessing  the  robustness  of  a  particular  control  law 
design  to  parameter  variations  and  unmodeled  dynamics  is  the  gain  margin.  This  is 
defined  as  the  extent  to  which  changes  can  be  inflicted  on  the  system  gain  without 
loss  in  stability,  to  this  end,  we  assume  that  the  heading  error  gain  K^  is  multiplied 
by  a  constant  C.  By  definition,  of  Hopf  bifurcation  occurs  when  a  pair  of  complex 
conjugate  eigenvalues  cross  into  the  right  hand  half-plane.  When  this  occurs  the 
system  will  deviate  from  a  steady  solution  in  an  oscillatory  manner.  This  deviation 
can  be  either  supercritical  or  subcritical  [Seydel  (1988)].  As  the  parameter  C  crosses 
the  critical  value,  one  pair  of  complex  conjugate  eigenvalues  of  the  linear  system 
matrix  crosses  transversely  the  imaginary  axis.  Locally,  as  C  approaches  Cent,  the 
periodic  solutions  are  located  on  the  two  dimensional  Eucledian  plane  spanned  by  the 
eigenvectors  of  the  Jacobian  matrix  of  the  system  which  corresponds  to  the  critical 
pair  of  eigenvalues.  In  this  chap  stability  properties  of  the  periodic  solutions  are 
established.  In  order  to  establish  those  properties  the  main  nonlinear  terms  that 
dominate  the  system  are  isolated.  Center  manifold  theory  is  used  to  reduce  the 
flow  to  a  two  dimensional  manifold.  The  method  of  averaging  is  then  applied  to  the 
reduced  system. 

The  critical  value  of  c  for  stability  of  straight  fine  motion  remains  the  same  as 
[Oral  (1993)],  which  is 

Ccrft  =  0.2658 


11 


This  is  because  the  dynamics  of  the  controller  are  independent  from  the  dynamics 
of  the  observer  as  explained  in  Chapter  II. 

B.     THIRD  ORDER  EXPANSIONS  OF  THE  SYSTEM  EQUATIONS 
1.     Perturbation  in  K<k 

In  the  previous  chapter  we  worked  on  the  linear  system.  Now  we  are  go- 
ing to  introduce  the  nonlinear  terms  to  our  compensator.  In  this  case  the  equations 
of  motion  are 

#    =    r  (3.1) 

r    =    ar  +  a3r3  +  b6  (3.2) 

y    =    sin*  (3.3) 

where 

6    =    <U'tanhf-^-J  (3.4) 

60    =    CK*(y-y)  +  Kr(r-r)  +  Kv{y-y)  (3.5) 

or  in  compact  form, 

X  =  f(x),         X=[*,  r,  y,  *,  f,  yf  (3.6) 

This  system  can  be  written  in  the  form 

X  =  AX  +  g{x)  (3.7) 

A  is  the  Jacobian  matrix  of  f(x)  evaluated  at  X  =  0,  and  g(x)  contains  all  nonlinear 
terms  of  Equations  (3.1),  (3.2),  and  (3.3).  Taylor  expansion  of  the  nonlinear  terms 
about  the  equilibrium,  where  we  keep  the  first  non-vanishing  nonlinear  coefficients 
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only,  gives 


simp     =    $__^3 
6 


6    =    60- 


1 


3£   ° 


(3.8) 
(3.9) 


'•at 


Substitution  of  Equations  (3.8)  and  (3.9)  into  Equations  (3.1),  (3.2),  and  (3.3)  gives 
us  the  A  matrix  in  Equation  (3.7)  as  follows, 


A  = 


0  10  0  0  0 

bcK*    a  +  bKr    bKy    -bcK*    -bKr    -bKy 


1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

— £q? 

0 

0 

0 

0 

a 

~ir 

0 

0 

0 

1 

0 

-4 

(3.10) 


The  nonlinear  parts  are, 


9{x)  = 


0 


a3r 


1  3SL° 

6 
0 

0 

0 


(3.11) 


If  we  introduce  the  transformation  matrix  (T)  of  eigenvectors  of  A  evaluated  at  the 
bifurcation  point, 


T    =    [rriij]        ij  =  1,2,3,4,5,6 


T'1    =    K]         i,j  =  1,2,3,4,5,6 


(3.12) 
(3.13) 


the  linear  change  of  coordinates, 


x  =  Tz ,         z  —  T    x 
transforms  system  (3.7)  into  its  normal  form 

z  =  T^ATz  +  T-xg(Tz) 


(3.14) 


(3.15) 
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With  this  transformation  we  get, 


T~XAT  = 


0 

—Wo 

0 

0 

0 

0 

w0 

0 

0 

0 

0 

0 

0 

0 

Pi 

0 

0 

0 

0 

0 

0 

Pi 

0 

0 

0 

0 

0 

0 

Pz 

0 

0 

0 

0 

0 

0 

Pa 

(3.16) 


Using  center  manifold  theory  we  get,  as  shown  in  [Chow  and  Mallet-Paret  (1977)], 


9  = 


0 

ill*}  +  ^22*?  +  ^23*?22  +  tuZiZ* 
*31*l  +  *»*!  +  *33*?*2  +  titZiZ* 

0 

0 

0 


Substitution  of  Equation  (3.14)  into  Equation  (3.7)  yields, 


(3.17) 


Z\    =    -w0z2  +  rnz\  +  rX2z\z2  +  T\*Z\z\  +  rXAz\ 


z2    =    w0zi  +r2\z\  +  r22z\z2  +  r^z^l  +  r2Az\ 


(3.18) 
(3.19) 


2.     Integral  Averaging 

We  write  Equations  (3.18)  and  (3.19)  in  the  form 


z\    =    -wQZ2  +  Fi(zuz2)i 

Z2      =     W0Zi+F2(zi1Z2) 


(3.20) 
(3.21) 


If  we  introduce  polar  coordinates  in  the  form, 


Z\  —  R  cos  0  ,         z2  =  R  sin  0 


(3.22) 


Equations  (3.20),  (3.21)  result  in 


R    =    F^R,  0)  cos  0  +  F2{R,0)  sin  0 
R0    =    wQR  +  F2(R,0)  cos  0-Fx{R,0)  sin  0 


(3.23) 
(3-24) 
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Equation  (3.23)  yields 

R  =  P{6)R?  (3.25) 

where  P(0)  is  a  27r-periodic  function  in  the  angular  coordinate  6.  If  Equation  (3.25) 
is  averaged  over  one  cycle  in  0,  we  get  an  equation  with  constant  coefficients, 

R  =  KR3  (3.26) 

where, 

K  =  ^-  r  P(0)  •  dO  (3.27) 

27T  JO 

Equation  (3.27)  is  simplified  after  evaluation  of  the  integral  as, 

K  =  g  [3rn  +  r13  +  r22  +  3r24]  (3.28) 

where  the  coefficients  are  as  follows, 

**11  =  7*12^21  +  ni3£3i 

r13  =  ^12^24  +  "13^34 

r22  =  n22^23  "H  n23^33 

r24  =  "22^22  +  "23^32 

where  the  ni2,  7113,  n22,  and  n^  are  the  elements  of  the  inverse  of  transformation 
matrix  T.  The  values  of  the  coefficients  £tJ  are  given  by  the  following  expressions 

l2\    =    a2m\l  -  be  [c3Klm34l  +  Klm\x  +  K Jm^  +  Z<?K%KTm2Almsl 

+3c  Kq,Kym4lm6i  +  3cAr#/^rm51m4i  +  3c/T#i^ym61m4i  -f  3KrKym61msi 
+ZK2Kym\xmQi  +  GcKy  KrKym4lm51m61 
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^22    =    a3m\2  -  bt  \czK%m\2  +  Kfm\2  +  ifjjm^  +  3c2  tfjirrm22  77153 

+3c2Kl Kym\2rn&2  +  3 cK* K?ml2m42  +  3cK^Klm\2mA2  -f  ZKTKlm\2m52 
Jr3KrKymb2me2  +  GcKy  KrKym42'^n^2rrH2 

£23    =    3a3m217n22  —  6/  3c  /Cj,m41m42  +  3^  m51m52  +  3/fy  m^  777,62 

+3c  KyjfKr  fm41 77152  +  277141 7714277151 J  +  3c  KyKy  (t7141  777-62  +  2771417714277161] 

+3cK<iiKl  (771^77142  +  2m5i m52m4i ]  +  3cKif,Ky  f m61m42  +  2m6im62m4iJ 

+3/^ Ky  ^77l61 77152  +  2m6i  77162^51  J  +  3KfKy  ^m51 77162  +  2771517715277261 J 
+6cKyKrKy  (771417715177162  +  (^4177152  +  "l42m5l)  m6l)J 

^24    =    3a3m21m22  —  6;  |3c  K^m^m^  +  3/{rm51m52  +  3Kym&xm62 

+3c  K\Kr  ^m42m5i  +  2m4im42m52]  +  3c  K^Ky  ^m42m6i  +  2m4im427n62j 
+ZcKyKl  (m22m4i  +  2m5im52m42]  +  ZcKyK*  (m62m4i  +  2m6im62m42] 
+ZKrKy  fm62m5i  +  2m6i 77162^52]  +  SKrKy  fm52m6i  +  2m5i 7775277162] 
+6cKq,KrKy  (771427715277161  +  (m4im52  +  m42m5i)  m62)J 


4i    =    _£mn 


^32  =  ~Fml2 

^33  =  -  «"hlm12 

^34  =  -2m"m12 

6,  =       b 


3*L 
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C.     RESULTS 

The  value  of  K  is  important  for  us  to  determine  whether  the  bifurcation  is 
supercritical  (K  <  0)  or  subcritical  (K  >  0).  In  this  study  we  wanted  to  see  the 
effects  of  observer  dynamics  to  our  system,  especially  the  value  of  K.  To  do  that 
we  used  the  Fortran  code  (Appendix  A)  for  the  numerical  results.  The  result  was 
that  the  value  of  K  was  not  affected  by  changes  in  the  observer  natural  frequency. 
The  reason  for  this  can  be  traced  back  to  the  definition  of  K,  Equation  (3.28).  It 
can  be  seen  that  in  the  definitions  for  r,j  and  l{j  only  the  first  two  eigenvectors  m,i, 
m,2  for  i  =  1,2,. .  .6  of  the  A  matrix,  Equation  (3.10)  appear.  Therefore,  we  have 
to  show  that  these  remain  the  same  for  all  observer  natural  frequencies. 

This  A  matrix,  Equation  (3.10),  actually  consists  of  4  block  matrices,  each 
3x3,  which  are  the  same  as  shown  in  Equation  (2.22).  Let  us  denote  the  A  matrix 
as  follows, 


A  = 


A    B 
0    C 


The  eigenvalues  of  A  can  be  computed  by 


A- XI        B 
0         C-XI 


=  0 


or 

\A  -  XI\  •  \C  -  XI\  =  0 

We  can  group  the  eigenvalues  in  two  different  groups:   Aitl-  for  i  =  1,2,3  are  the 
eigenvalues  of  A  an  A2l,  for  i  =  1,2,3  are  the  eigenvalues  of  C.    The  eigenvectors 
associated  with  these  eigenvalues  can  be  found  as  follows. 
For  A  =  Alt,, 


A  -  XUI         B 

Vl 

0 

0          C  -  XUI 

v2 

0 
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and 

[«4-Ai,,-/]H  +  [3][t;2]    =    0 

[0]N  +  [C-A1)t/][t;2]     =    0 

Since  A1(,  is  an  eigenvalue  of  A  and  the  eigenvalues  of  A  and  C  are  distinct,  the 
matrix  [C  —  \\,%I\  is  nonsingular  which  means  that  fa)  =  0.  Therefore,  we  get 

[^-a1,,/]M  =  o 

which  means  that  V\  is  an  eigenvector  of  A.  Therefore,  the  first  three  eigenvectors  of 
A  are  the  same  as  the  eigenvectors  of  A  and  they  are  independent  of  the  dynamics 
of  the  observer.  Of  course,  the  remaining  three  eigenvectors  of  A  depend  on  the 
oberver  natural  frequency,  but,  as  we  pointed  out  earlier,  none  of  these  appear  in 
the  definition  of  the  nonlinear  stability  coefficient  K. 
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IV.  COMPENSATOR  IN  A  DIFFERENT  BASIS 

A.  CRITICAL  VALUE  OF  C 

If  we  look  at  Equation  (3.6),  we  can  see  that  the  basis  for  our  system  was 
[X,  X].  Now  we  are  going  to  represent  our  system  in  [X,X]  basis  where  X  is 
the  estimate  of  X.  In  this  compensator  basis  the  critical  value  of  C  in  Equation 
(4.3)  is  no  longer  constant.  Therefore,  we  used  a  Fortran  code  (Appendix  B)  to 
calculate  the  critical  C  values  for  different  observer  natural  frequencies.  A  plot 
of  these  critical  C  values  for  different  observer  natural  frequencies  can  be  seen  in 
Fig.  4.1.  The  observer  natural  frequencies  are  in  nondimensional  seconds  whereas 
the  control  natural  frequencies  are  normalized  with  respect  to  the  corresponding 
observer  natural  frequencies.  The  system  is  unstable  for  values  of  C  less  than  the 
critical  value. 

B.  THIRD  ORDER  EXPANSIONS  OF  THE  SYSTEM  EQUATIONS 

1.     Perturbation  in  Kq, 

In  the  previous  chapters  we  worked  on  the  [X,  X]  basis.  Now  we  are 
going  to  represent  our  system  in  the  new  basis,  which  is  [X,  X],  where  X  is  the 
estimate  of  X.  Our  equations  of  motion  were, 

i    =    r 

r    =    ar  +  a^r3  +  b6  (4.1) 

y    =    sin  $ 
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CCRTT  vs   WN(normaiized) 


0.35 


0.15- 


1.6  1.8  2  12  2.4  16  2.3  3 

WN(normaiized) 


Figure  4.1:     C^  versus  natural  frequency  for  K^. 


where 


5    =    6«At  tanh  I  - — 


J*a£ , 


So    =     Ctf*#  +  /Cr  +  jrvy 


(4.2) 
(4.3) 


or  in  compact  form, 


X  =  f(x),  X  =    *,r,y,*,r,y 


This  system  can  be  written  in  the  form 


(4.4) 


X  =  AX  +  g{x) 
20 


(4.5) 


A  is  the  Jacobian  matrix  of  f(x)  evaluated  at  X  =  0,  and  g(x)  contains  all  non- 
linear terms  of  Equation  (4.1).  Taylor  expansion  of  the  nonlinear  terms  about  the 
equilibrium,  where  we  keep  the  first  non-vanishing  nonlinear  coefficients  only,  gives 

(4.6) 

(4.7) 

Substitution  of  Equations  (4.6)  and  (4.7)  into  Equation  (4.1)  gives  us  the  A  matrix 
in  Equation  (4.5)  as  follows 

0 


sin 

tf 

= 

*- 

i*3 
6 

6 

= 

*o- 

1 

A  = 


0    10        0         0 

0    a     0  bcK*  bKr  bKy 

10     0        0         0  0 

0    0    l«       0          1  -l* 

0    0£r  bcK*  bKr  -tr  +  bKy 


0    0£, 


1 


0 


-L 


(4.8) 


The  nonlinear  parts  are, 


0(*)  = 


a3r3  - 


ZSL 


-'sat 

_i#3 

6 
0 

4  * 


ML   ° 


JSAt 

0 


(4.9) 


If  we  introduce  the  transformation  matrix  (T)  of  eigenvectors  of  A  evaluated  at  the 
bifurcation  point, 


T    =    [rriij]        i,;  =  1,2,3,4,5,6 
T-1    =    K-]        i,j  =  1,2,3,4,5,6 


(4.10) 
(4.11) 


the  linear  change  of  coordinates, 


x  =  Tz  ,         z  —  T    x 


(4.12) 
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transforms  system  (4.5)  into  its  normal  form 

z  =  T~xATz  +  T'lg(  Tz) 
At  the  bifurcation  point 


(4.13) 


T~lAT  = 


(4.14) 


0  -w0  0  0  0  0 

w0      0  0  0  0  0 

0       0  P1  0  0  0 

0       0  0  P2  0  0 

0       0  0  0  P3  0 

0       0  0  0  0  P4  . 

with  wo  >  0  and  P,  <  0. 

Using  center  manifold  theory  we  get,  as  shown  in  [Chow  and  Mallet-Paret 

(1977)], 

0 
'21*?  +  '22*?  +  t2Zz\z2  +  l2Azxz\ 

4l*i  +  '32*2  +  '33*1*2  +  '34*1*2 
0 

'51*1  +  '52*2  +  '53*i*2  +  '54*1*2 

0 

Substitution  of  Equation  (4.12)  into  Equation  (4.5)  yields, 


9  = 


(4.15) 


ii     =    -w0z2  +  rnzf  +  ri2zjz2  +  ri3*i*f  +  f*i4*2 


i2    =    w0zi  +  r2i  2?  +  r22Zi2*2  +  ^23*1  z\  +  r2Az\ 


(4.16) 
(4.17) 


where  the  terms  rfJ  are  evaluated  by  a  Fortran  code  (Appendix  C). 

For  values  of  C  close  to  its  critical  value,  Equations  (4.16)  and  (4.17) 
become, 

ii    =    a'ezx  -  (wo  +  w'e)z2  +  rnz\  +  r^zf z2  +  ri32i2|  +  r^zf        (4.18) 
z2    =    («>o  +  ^'e)zi  +<*'£z2  +  r2lz\  +  r22z\z2  +  r23zxz\  +  r2Az\        (4.19) 

where  e  is  the  difference  between  C  and  its  critical  value  C*.  The  terms  a'  and  u/ 
denote  the  derivative  of  the  real  and  imaginary  part,  respectively,  of  the  critical  pair 
of  the  eigenvalues  with  respect  to  C  evaluated  at  C* 
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2.     Integral  Averaging 

We  write  Equations  (4.18)  and  (4.19)  in  the  form, 

ix    =    a'ezx  -  (w0  +  w'e)z2  +  Fx(zi,z2)  (4.20) 

z2    =    (u>o  +  w'e)z1  +  a'ez2  +  i^i,  22)  (4.21) 

If  we  introduce  polar  coordinates  in  the  form, 

zl  =  Rcos0,         z2  =  Rs'm$  (4.22) 

Equations  (4.20),  (4.21)  result  in 

R    =    a'eR  +  Fx  {R,  9)  cos  0  +  F2(R,0)  sin  0  (4.23) 

R6    =    (w0  We^  +  i^W^cosfl-F^i^sinfl  (4.24) 

Equation  (4.23)  yields 

R  =  a'eR  +  P^)/*3  (4.25) 

where  P(0)  is  a  2x-periodic  function  in  the  angular  coordinate  6.  If  Equation  (4.25) 
is  averaged  over  one  cycle  in  $  [Chow  and  Mallet-Paret  (1977)],  we  get  an  equation 
with  constant  coefficients, 

R  =  a'eR  +  KR3  (4.26) 

where 

K  =  ±J*P(e)dO  (4.27) 

Equation  (4.27)  is  simplified  after  evaluation  of  the  integral  as, 

K  =  g  Pfu  +  r13  +  r22  +  3r24]  (4.28) 

where  the  coefficients  are  as  follows 
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1*13  =  "12^24  +  "13^34  +  Kl^SA 
r22  =  ^22^23  +  7*23^33  +  **25^53 
1"24      =      "22^22  +  "23^32  +  "25^52 

where  the  n^,  7213,  72i5,  7222,  «23»  and  7225  are  the  elements  of  the  inverse  of  trans- 
formation matrix  T.  The  values  of  the  coefficients  £21,  £22>  ^23>  4*4>  ^31  ?  ^32?  ^33>  and 
£34  are  the  same  as  in  Chapter  III.  The  values  of  the  other  £{j  coefficients  are  given 
by  the  following  expressions: 

hi    =    -be  [c3I<lm341  +  Kzm\x  +  ifjmj,  +  3c2KlKrm241msl 

+Zc2  K%Kym24lm6i  +  ZcK^K^Tn^m^  +  "icK>iK2m\xmAx  +  S/^r^mgims! 

+3A'r.ftrym517726i  4-  ^cK^KTKymA\m^\ms\\ 

^52    =    —  bt  \c  K^m42  +  KrmS2  +  ^ym62  "I"  3c  KiffKrm42ms2 

+3c  KqKym42mQ2  +  3cKyKrm52m42  +  ScA^A^  772^77242  +  3KrKym62ms2 

+3KrKym52m62  +  6c/^/C^y*™42m52rr*62j 


^sa    =    —6, 


3c  A^m,,! 777,42  +  3Arrm51m52  +  3A"ym61m62 
+3c  KyKr  [rn4lms2  +  2m41m42m51J  +  3c  K9Ky  fm41m62  +  2m41m42m6iJ 

+3cKyK?  (7725,77142  +  27715x771527724!)  +  3cKyK*  (jn^m^  +  277261m6277241  J 
+3/C^y  ("^6lm52  +  2m61 772627725! J  +  3K?Ky  (rn\xm62  +  2772517725277261J 
+6cKif,KrKy  (77241 77251 77262  +  ("^41^52  +  "^TTlsi)  7726l )] 
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£54    =    —  be  1 3c  K^rri4\mA2  +  tiKTms\mb2  +  3Arym6im62 

+$<?K%Kr  (mj2m5i  +  2m41m42m52)  +  Z<?K%Ky  (rn\2rn^i  +  2m4im42m62J 
+3cKq,K?  (rn\2mAi  +  2m5im52m42)  +  3cK*Ky  (m\2mAi  +  2m6im62ra42J 
+3/fr#y  (m62m51  +  2m61m62m52)  +  3K?Ky  (rn\2m6l  +  2m51m52m62J 
+6cif*  A^-ft^  (m42m52m6i  +  (m4im52  +  m42m51)  m^)] 

C.     RESULTS 

Existence  and  stability  of  limit  cycles  can  be  determined  by  analyzing  the 
equilibrium  points  of  the  averaged  Equation  (4.26),  which  correspond  to  periodic 
solutions  in  Zi,  z2  as  can  be  seen  from  Equation  (4.22).  From  Equation  (4.26)  we 
can  easily  see  that: 

1.  If  a'  >  0,  then 

(a)  if  K  >  0,  then  unstable  periodic  solutions  co-exist  with  the  stable  equi- 
librium for  e  <  0,  and 

(b)  if  K  <  0,  then  stable  periodic  solutions  co-exist  with  the  unstable  equi- 
librium for  e  >  0. 

2.  If  a'  <  0,  then 

(a)  if  K  >  0,  then  unstable  periodic  solutions  co-exist  with  the  stable  equi- 
librium for  £  >  0,  and 

(b)  if  K  <  0,  then  stable  periodic  solutions  co-exist  with  the  unstable  equi- 
librium for  £  <  0. 
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We  refer  to  K  <  0  as  the  supercritical,  and  K  >  0  as  the  subcritical  PAH  bi- 
furcation. In  the  supercritical  case,  after  the  equilibrium  state  loses  its  stability 
the  system  converges  to  a  stable  periodic  solution  with  amplitude  which  increases 
continuously  as  the  difference  e  is  increased. 

In  the  subcritical  case,  however,  before  the  equilibrium  state  loses  stability,  its 
domain  of  attraction  becomes  very  small  since  it  is  bounded  by  the  amplitudes  of  the 
unstable  limit  cycles.  In  such  a  case,  an  initial  disturbance  of  sufficient  magnitude 
can  throw  the  system  off  the  nominal  path  even  before  its  domain  of  attraction 
has  completely  shrunk  to  zero.  As  the  nominal  equilibrium  becomes  unstable,  the 
system  jumps  to  a  different  state  of  motion  with  a  locally,  at  e  =  0,  discontinuous 
increase  in  the  amplitude  [Papoulias  (1993)]. 

In  our  case,  the  value  of  a'  is  always  negative,  which  can  be  seen  easily  from 
Figure  4.1.  If  we  look  at  the  nature  of  the  curve  for  the  critical  value  of  C  for 
different  natural  frequencies,  we  will  see  that  as  the  value  of  critical  C  decreases  for 
the  same  natural  frequency,  the  system  becomes  unstable. 

After  using  a  Fortran  code  (Appendix  C)  we  observed  that  the  nonlinear  sta- 
bility coefficient  K  depends  on  observer  dynamics.  Figure  4.2  shows  that  for  a  given 
control  design,  the  observer  must  be  as  responsive  as  possible  to  ensure  negative  K 
(stable  limit  cycle).  On  the  other  hand,  for  a  given  observer  design,  if  the  control 
law  is  too  slow  we  get  subcritical  behavior  (unstable  limit  cycle). 
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Figure  4.2:     Kr    versus  wn  for  different  observer  wn. 


27 


[THIS  PAGE  INTENTIONALLY  LEFT  BLANK] 


28 


V.  CONCLUSIONS  AND  RECOMMENDATIONS 

A.     CONCLUSIONS 

An  investigation  of  the  nonlinear  dynamic  response  characteristics  of  a  marine 
vehicle  has  been  presented.  Particular  emphasis  in  this  work  was  placed  on  analyzing 
the  effects  of  observer  design  on  system  response  after  initial  loss  of  stability  of 
straight  line  motion.  Bifurcation  theory  techniques  were  utilized  in  order  to  assess 
that  behavior.  The  main  conclusions  of  this  work  can  be  summarized  as  follows. 

1.  There  exists  a  critical  point  for  a  certain  combination  of  system  gains  and 
system  parameters  for  stability  of  straight  line  motion.  The  loss  of  stability 
occurs  generically  in  the  form  of  Poincare-Andronov-Hopf  bifurcations.  As 
the  parameter  crosses  its  critical  value,  a  family  of  periodic  orbits,  self  sus- 
tained oscillations  develops.  Center  manifold  reduction  and  integral  averaging 
techniques  were  used  in  order  to  establish  the  direction  of  the  bifurcation  and 
stability  of  the  resulting  periodic  solutions  [Papoulias,  Oral  (1993)]. 

2.  For  [X,  X]  basis  the  critical  point  does  not  depend  on  the  observer  dynamics 
(separation  principle).  The  nonlinear  stability  coefficient  K  was  not  influenced 
by  observer  dynamics.  The  previous  reduction  process  shows  that  K  depends 
on  the  first  two  eigenvectors  of  the  6x6  matrix  A.  Matrix  algebra  shows  that 
these  eigenvectors  are  associated  only  with  the  controller  dynamics. 

3.  For  [X,  X)  basis  the  critical  value  depends  on  observer  dynamics.  For  a  given 
control  design,  the  observer  must  be  as  responsive  as  possible  to  maximize  the 
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region  of  stability.  On  the  other  hand,  for  a  given  observer  design,  the  control 
must  be  as  slow  as  possible  to  maximize  region  of  stability. 

4.  The  nonlinear  stability  coefficient  K  depends  on  observer  dynamics  for  this 
basis.  For  a  given  control  design,  the  observer  must  be  as  responsive  as  possible 
to  ensure  negative  K  (stable  limit  cycles).  In  this  benign  loss  of  stability 
the  resulting  periodic  solutions  are  continuous  single-valued  functions  of  the 
parameter  distance  from  its  critical  value.  On  the  other  hand,  for  a  given 
observer  design,  if  the  control  law  is  too  slow  we  get  sub  critical  behavior 
(unstable  limit  cycles).  In  such  a  case,  the  periodic  solutions  develop  with 
what  appears  to  be  a  discontinuous  increase  in  the  amplitude  of  oscillations 
[Papoulias,  Oral  (1993)]. 

B.     RECOMMENDATIONS 

The  differences  between  the  two  bases  with  respect  to  robustness  properties  of 
the  system  have  to  be  analyzed. 
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APPENDIX  A 
HOPF  BIFURCATION  PROGRAM  FOR  [X,X]  BASIS 


PROGRAM  HTKPSI 
C      HOPF  BIFURCATIONS 
C      NOMOTO'S  FIRST  ORDER  MODEL 
C 

C      CALCULATIONS  FOR  K  AND  CCRITICAL  IF  KPSI  CHANGES  W/  C 
C 
C234567 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
C 

DOUBLE  PRECISION  Kl , K2 , K, LPSI , LY, LR, IZ, L, 
&  MASS ,  NV, NR, NVDOT, NRDOT,  NDRS ,  NDRB,  KPSI ,  KR,  KY,  K3  , 

Ec  L21  ^22^23^24^31^32^33^34^51^52^53^54, 

&  Ml  1 ,  Ml  2  ,  Ml  3  ,  Ml  4 ,  Ml  5 ,  Ml  6 ,  M2 1 ,  M2  2 ,  M2  3 ,  M2  4 ,  M2  5 ,  M2  6 , 

Sc  M31 ,  M32  ,  M33  ,  M34 ,  M3 5 ,  M3 6 ,  M41 ,  M42 ,  M43  ,  M44  ,  M45 , M46 , 

&  M51 , M52  ,  M53  , M54  ,  M55 ,  M56 ,  M61 ,  M62 ,  M63  ,  M64 ,  M65 ,  M66 , 

Sc  Nil ,  N12  ,  N13  ,  N14 ,  N15 ,  N16 ,  N2 1 ,  N22 ,  N23  ,  N24 ,  N2 5 ,  N2 6 , 

&  N3 1 ,  N3  2  ,  N3  3  ,  N3  4 ,  N3  5 ,  N3  6 ,  N4 1 ,  N4  2 ,  N4  3 ,  N4  4 ,  N4  5 ,  N4  6 , 

Sc  N51 ,  N52 ,  N53  ,  N54 ,  N55 ,  N56 ,  N61 ,  N62  ,  N63  ,  N64 ,  N65 ,  N66 

C 

DIMENSION  AMAT (6,6),T(6,6), TINV (6,6), FV1 ( 6 ) , IV1 ( 6 ) , YYY (6,6) 

DIMENSION  WR  (  6  )  ,  WI  (  6 )  ,  TSAVE  (6,6),  TLUD  (6,6),  IVLUD  (  6  )  ,  SVLUD  (  6 ) 

DIMENSION  ASAVE ( 6 , 6 ) , Al ( 6 , 6 ) , A2 ( 6 , 6 ) 

OPEN    (11,FILE='AKPSI.MAT' ,STATUS=/NEW ) 
C 

WEIGHT=435.0 

IZ     =45.0 

L      =7.3 

RHO    =1 . 94 

G      =32.2 

XG     =0.0104 

MASS   =WEIGHT/G 

MASS   =MASS/ (0.5*RHO*L**3) 

IZ     =IZ/ (0.5*RHO*L**5) 

XG     =XG/L 

YRDOT  =-0.00000 

YVDOT  =-0.03430 

YR     =+0.00000 

YV     =-0.10700 

YDRS   =+0.01241 

YDRB   =+0.01241 

NRDOT  =-0.00047 

NVDOT  =-0.00000 

NR     =-0.00390 
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NV     =-0.00000 

NDRS   =-0.337*YDRS 

NDRB   =+0.283*YDRS 

DH   = (IZ-NRDOT) * (MASS-YVDOT) - 
Sc  (MASS*XG-YRDOT)  *  (MASS*XG-NVDOT) 

All= ( ( IZ-NRDOT) *YV- (MASS*XG-YRDOT) *NV) /DH 

A12=( (IZ-NRDOT) * (-MASS+YR) - 
8c  (MASS*XG-YRDOT)  *  (  -MASS*XG+NR)  )  /DH 

A21= ( (MASS-YVDOT) *NV- (MASS*XG-NVDOT) *YV) /DH 

A22=( (MASS-YVDOT) * ( -MASS*XG+NR) - 
&      (MASS*XG-NVDOT) * ( -MASS+YR) ) /DH 

Bll=( (IZ-NRDOT) *YDRS-(MASS*XG-YRDOT) *NDRS) /DH 

B12  = ( (IZ-NRDOT) *YDRB- (MASS*XG-YRDOT) *NDRB) /DH 

B21= ( (MASS-YVDOT) *NDRS- (MASS*XG-NVDOT) *YDRS) /DH 

B22= ( (MASS-YVDOT) *NDRB- (MASS*XG-NVDOT) *YDRB) /DH 

Bl  =B11-B12 

B2  =B21-B22 
C 

200  WRITE  (*,1004) 

READ   (*,*)     WNMIN/WNMAX,IWN 

INCR=IWN 

WRITE  ( * , * )  ' ENTER  OBSERVER  WN ' 

READ  [*,*)    WNO 
205  WRITE(*,1007) 

READ  ( * , * )    A3 
C     DO  is  Dsat 
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WRITE 

(*,1006) 

READ 

(*,*)     DO 

WRITE 

(*,1008) 

READ 

[*,*)       CCRIT 

c 
c 


Cl= (A11*A22-A21*A12) * (A21*B1-A11*B2 ) 

C2=(A11+A22) * (A21*B1-A11*B2)+B2* (A11*A22-A21*A12 ) 

C3=- (A21*B1-A11*B2) **2 

A=C1/C2 

B=C3/C2 

DO  1  11=1, INCR 


WN   =WNMIN+ (WNMAX-WNMIN) * (II-l) / (INCR-1 

print  *,wn 
ALPHA0=WN**3 
ALPHA1=2.15*WN**2 
ALPHA2=1.75*WN 
KPSI=-ALPHA1/B 
KY   =-ALPHA0/B 
KR   =- (ALPHA2+A) /B 
GAMA0=WNO**3 
GAMA1=2 .15*WNO**2 
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GAMA2  =  : 

L.75*WNO 

LY=A+GAMA2 

LPSI=A*LY+GAMA1 

LR=A*LPSI+GAMAO 

C234567 

C           A           MATRIX 

AMAT(1,1) =0.0 

AMATd 

2)=1.0 

AMATd 

3)=0.0 

AMATd 

4)=0.0 

AMATd 

5)=0.0 

AMATd 

6)=0.0 

AMAT(2 

, 1)=B*CCRIT*KPSI 

AMAT(2 

2)=A+(B*KR) 

AMAT(2 

3) =B*KY 

AMAT(2 

, 4)=-B*CCRIT*KPSI 

AMAT(2 

5)=-B*KR 

AMAT(2 

6) =-B*KY 

AMAT(3 

1)=1.0 

AMAT(3 

2)=0.0 

AMAT(3 

3)=0.0 

AMAT(3 

4)=0.0 

AMAT(3 

5)=0.0 

AMAT(3 

,6)=0.0 

AMAT(4 

,1)=0.0 

AMAT(4 

,2)=0.0 

AMAT(4 

3)=0.0 

AMAT(4 

,4)=0.0 

AMAT(4 

5)=1.0 

AMAT(4 

6)=-LPSI 

AMAT(5i 

1)=0.0 

AMAT(51 

2)=0.0 

AMAT(5 

3)=0.0 

AMAT(5 

4)=0.0 

ZMAT{5, 

5)=A 

AMAT(5 

6) =-LR 

AMAT(6, 

1)=0.0 

AMAT(6, 

2)=0.0 

AMAT(6 

3)=0.0 

AMAT(6 

4)=1.0 

AMAT(6 

5)=0.0 

AMAT(6 

6) =-LY 

DO    11    1=1, 6 

DO    12    J=l,6 

ASAVE ( I , J) =AMAT ( I , J) 

12                 CONTINUE 

11           CONTINUE 

CALL   RG(6, 6,AMAT,WR,WI , 1,YYY, IV1 , FV1 , 

I  ERR) 

CALL    DSOMEG(IEV,WR,WI, OMEGA, CHECK) 

C                WRITE 

(*,*)     IEV 
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WRITE  (60,*)   (WR(IREAL) ,IREAL=1, 6) 

OMEGA0=OMEGA 

DO  5  1=1, 6 

T(I, 1) =YYY(I, IEV) 

T(I,2) =-YYY(I, IEV+1) 
5      CONTINUE 

IF(IEV.EQ.l)  GO  TO  13 
IF(IEV.EQ.2)  GO  TO  14 
IF(IEV.EQ.3)  GO  TO  15 
IF(IEV.EQ.4)  GO  TO  16 
IF(IEV.EQ.5)  GO  TO  17 
STOP  3004 

13  DO  21  1=1,6 

T ( I , 3 ) =YYY (1,3) 
T ( I , 4 ) =YYY (1,4) 
T ( I , 5 ) =YYY (1,5) 
T  ( 1 ,  6  ) =YYY (1,6) 

2 1  CONTINUE 
GO  TO  3  0 

14  DO  22  1=1, 6 

T ( I , 3 ) =YYY (1,1) 
T ( I , 4 ) =YYY (1,4) 
T ( I , 5 ) =YYY (1,5) 
T ( I , 6 ) =YYY (1,6) 

22  CONTINUE 
GO  TO  3  0 

15  DO  23  1=1,6 

T ( I , 3 ) =YYY (1,1) 
T ( I , 4 ) =YYY (1,2) 
T ( I , 5 ) =YYY (1,5) 
T ( I , 6 ) =YYY (1,6) 

23  CONTINUE 
GO  TO  3  0 

16  DO  24  1=1, 6 

T ( I , 3 ) =YYY (1,1) 
T ( I , 4 ) =YYY (1,2) 
T ( I , 5 ) =YYY (1,3) 
T ( I , 6 ) =YYY (1,6) 

24  CONTINUE 
GO  TO  3  0 

17  DO  25  1=1,6 

T ( I , 3 ) =YYY (1,1) 
T ( I , 4 ) =YYY (1,2) 
T ( I , 5 ) =YYY (1,3) 
T ( I , 6 ) =YYY (1,4) 

2  5     CONTINUE 

3  0     CONTINUE 
C 

C      NORMALIZATION  OF  THE  CRITICAL  EIGENVECTOR 
C 
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CALL  NORMAL (T) 

INVERT  TRANSFORMATION  MATRIX 

DO  2  1=1, 6 
DO  3  J=l, 6 

TINV(I, J) =0.0 
TSAVE(I, J)=T(I, J) 
CONTINUE 
CONTINUE 

CALL  DLUD (6,6, TS AVE , 6 , TLUD , IVLUD ) 
DO  4  1=1, 6 

IF  (IVLUD(I) .EQ.O)  STOP  3003 
CONTINUE 

CALL  DILU (6,6, TLUD , IVLUD , SVLUD ) 
DO  8  1=1, 6 
DO  9  J=l, 6 

TINV(I, J)=TLUD(I, J) 
CONTINUE 
CONTINUE 

CHECK   Inv(T)*A*T 

IMULT=1 

IF  (IMULT.EQ.l)  CALL  MULT (TINV, ASAVE, T, A2 ) 

IF  (IMULT.EQ.O)  STOP  3  007 

P1=A2 (1,1) 

P2=A2 (2,2) 

P=A2(3,3) 

Q=A2(4,4) 

R=A2 (5,5) 

S=A2 (6,6) 

WRITE (21, *)P1,P2,P,Q/R,S 

DEFINITION  OF  Nij 

N11=TINV(1,1) 
N21=TINV(2,1) 
N31=TINV(3,1) 
N41=TINV(4,1) 
N51=TINV(5,1) 
N61=TINV(6,1) 
N12=TINV(1,2) 
N22=TINV(2,2) 
N32=TINV(3,2) 
N42=TINV(4,2) 
N52=TINV(5,2) 
N62=TINV(6,2) 
N13=TINV(1,3) 
N23=TINV(2,3) 
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N33=TINV(3,3) 
N43=TINV(4,3) 
N53=TINV(5,3) 
N63=TINV(6,3) 
N14=TINV(1,4) 
N24=TINV(2,4) 
N34=TINV(3,4) 
N44=TINV(4,4) 
N54=TINV(5,4) 
N64=TINV(6,4) 
N15=TINV(1,5) 
N25=TINV(2,5) 
N3  5=TINV(3,5) 
N45=TINV(4/5) 
N55=TINV(5, 5) 
N65=TINV(6,5) 
N16=TINV(1, 6) 
N2  6=TINV(2, 6) 
N3  6=TINV(3,6) 
N46=TINV(4, 6) 
N56=TINV(5,6) 
N66=TINV(6, 6) 

C 

C  DEFINITION   OF   Mij 

C 

M11=T(1,1) 
M21=T(2,1) 
M31=T(3/1) 
M41=T(4,1) 
M51=T(5,1) 
M61=T(6,1) 
M12=T(1,2) 
M22=T(2,2) 
M32=T(3,2) 
M42=T(4,2) 
M52=T(5,2) 
M62=T(6,2) 
M13=T(1,3) 
M23=T(2,3) 
M33=T(3,3) 
M43=T(4,3) 
M53=T(5,3) 
M63=T(6,3) 
M14=T(1,4) 
M24=T(2,4) 
M34=T(3,4) 
M44=T(4,4) 
M54=T(5,4) 
M64=T(6,4) 
M15=T(l/5) 
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M25=T(2,5) 

M35=T(3,5) 

M45=T(4,5) 

M55=T(5,5) 

M65=T(6,5) 

M16=T(1,6) 

M2  6=T(2, 6) 

M3  6=T(3,6) 

M46=T(4, 6) 

M56=T(5,6) 

M66=T(6,6) 

WRITE (70, *) 

N11,N12,N13 

WRITE (71, *) 

N14,N15,N16 

WRITE (72, *) 

N21,N22,N23 

WRITE (73,*) 

N24,N25,N2  6 

WRITE (74, *) 

N31,N32,N33 

WRITE (75,*) 

N34,N35,N36 

WRITE (76, *) 

N41,N42,N43 

WRITE (77, *) 

N44,N45,N46 

WRITE(78, *) 

N51,N52,N53 

WRITE (79, *) 

N54,N55,N56 

WRITE (80, *) 

N61,N62,N63 

WRITE (81, *) 

N64,N65,N66 

C       Kl=l./8.* ( (ALPHA2**3)+ALPHA0) / (ALPHA2) 

C       K2=3.*A3-. 

5* (ALPHA2**2) /ALPHAO 

C       K3=l./ ( (B**2) *(D0**2) ) * (ALPHA2+A) * ( (ALPHA0/ALPHA2 ) + ( A**2 

C       K=K1*(K2+K 

C 

C      print  * ,  wj 
p 

3) 

i, k, ccrit 

BL=B/ (3*D0**2) 
C 
C2345 6789012345678901234567890123456789 012345678901234567890123 

6789012 

L21=A3*M21**3-BL* ( CCRIT* *3*KPSI**3*M41**3+KR**3*M51**3+ 
&  KY**3*M61**3+ 

&  3*CCRIT**2*KPSI**2*KR*M41**2*M51+ 

&  3*CCRIT**2*KPSI**2*KY*M41**2*M61+ 

&  3*CCRIT*KPSI*KR**2*M51**2*M41+3*CCRIT*KPSI*KY**2*M61**2*M41+ 
&  3*KR*KY**2*M61**2*M51+3*KR**2*KY*M51**2*M61+ 

&  6*CCRIT*KPSI*KR*KY*M41*M51*M61) 

L22=A3*M22**3-BL* ( CCRIT* *3*KPSI**3*M42**3+KR**3*M52**3+ 
Sc  KY**3*M62**3  + 

&  3*CCRIT**2*KPSI**2*KR*M42**2*M52+ 

&  3*CCRIT**2*KPSI**2*KY*M42**2*M62+ 

&  3*CCRIT*KPSI*KR**2*M52**2*M42+3*CCRIT*KPSI*KY**2*M62**2*M42+ 
Sc  3*KR*KY**2*M62**2*M52  +  3*KR**2*KY*M52**2*M62  + 

Sc  6*CCRIT*KPSI*KR*KY*M42*M52*M62) 

L2  3  =3  * A3  *M2 1**2  *M2 2 -BL* ( 3  *CCRIT*  * 3  *KPSI *  *  3  *M4 1**2  *M42  + 
Sc  3*KR**3*M51**2*M52  + 
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&  3*KY**3*M61**2*M62+ 

&  3*CCRIT**2*KPSI**2*KR* (M41**2*M52+2*M41*M42*M51 ) + 

&  3*CCRIT**2*KPSI**2*KY* (M41**2 *M62+2*M41*M42*M61 ) + 

Sc  3*CCRIT*KPSI*KR**2*  (M51**2*M42+2*M51*M52*M41 )  + 

Sc  3*CCRIT*KPSI*KY**2*  (M61**2*M42+2*M61*M62*M41 ) + 

&  3*KR*KY**2*  (M61**2*M52  +  2*M61*M62*M51 )  + 

Sc  3*KR**2*KY*  (M51**2*M62  +  2*M51*M52*M61)  + 

Sc         6*CCRIT*KPSI**KR*KY*  (M41*M51*M62+  (M41*M52+M42*M51)  *M61)  ) 

L24=3*A3*M21*M22**2-BL* (3*CCRIT**3*KPSI**3*M41*M42**2+ 
Sc  3*KR**3*M51*M52**2  + 

Sc  3*KY**3*M61*M62**2  + 

Sc  3*CCRIT**2*KPSI**2*KR*  (M42**2*M51+2*M41*M42*M52  )  + 

Sc  3*CCRIT**2*KPSI**2*KY*  (M42**2*M61  +  2*M41*M42*M62  )  + 

Sc  3*CCRIT*KPSI*KR**2*  (M52**2*M41  +  2*M51*M52*M42  )  + 

&  3*CCRIT*KPSI*KY**2* (M62**2*M41+2*M61*M62*M42 ) + 

Sc  3*KR*KY**2*  (M62**2*M51+2*M61*M62*M52)  + 

&  3*KR**2*KY* (M52**2*M61+2*M51*M52*M62) + 

&  6*CCRIT*KPSI*KR*KY* (M42*M52*M61+ (M41*M52+M42*M51 ) *M62 ) ) 

L31=(-l./6. )*M11**3 

L32=(-l./6. )*M12**3 

L33=(-l./2. )*M11**2*M12 

L34=(-l./2. ) *M11*M12**2 

R11=N12*L21+N13*L31 

R12=N12*L23+N13*L33 

R13=N12*L24+N13*L34 

R14=N12*L22+N13*L32 

R21=N22*L21+N23*L31 

R22=N22*L23+N23*L33 

R23=N22*L24+N23*L34 

R24=N22*L22+N23*L32 

K=(3*R11+R13+R22+3*R24) /8 
WRITE  (11,2001)  WN,K,CCRIT 
.  CONTINUE 


STOP 

1004 

FORMAT 

('  ENTER  MIN,MAX,  AND  INCREMENTS  OF 
WN(stepstodivide) ' ) 

1006 

FORMAT 

( '  ENTER  DELTAS AT . ' ) 

1007 

FORMAT 

( '  ENTER  A3 ' ) 

1008 

FORMAT 

( '  ENTER  CCRIT' ) 

2001 

FORMAT 
END 

(3E15.5) 

SUBROUTINE  DSOMEG ( UK, WR, WI , OMEGA, CHECK) 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
DIMENSION  WR ( 6 ) , WI ( 6 ) 
CHECK=-1.0D+25 
DO  1  1=1, 6 

IF  (WR(I) .LT. CHECK)  GO  TO  1 

CHECK=WR(I) 
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IJ=I 

CONTINUE 

OMEGA=DABS (WI ( I J) ) 

IF  (WI(IJ) .GT.O.DO)  IJK=IJ 


IF  (WI(IJ) .LT.O.DO)  IJK=IJ-1 

RETURN 

END 


SUBROUTINE  NORMAL (T) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  T ( 6 , 6) , TNOR (6,6) 

CFAC=T(1,1) **2+T(l,2) **2 

IF  (DABS(CFAC) . LE . (l.D-10) )  STOP  4001 

TNOR (1,1)=1. DO 

TNOR(2,l)=(T(l,l) *T(2,1)+T(2,2) *T(1,2)  ) /CFAC 

TNOR (3,1)=(T(1,1)*T(3,1)+T(3,2)*T(1,2)) /CFAC 

TNOR (4,1)=(T(1,1)*T(4,1)+T(4,2)*T(1,2)) /CFAC 

TNOR (5,1)=(T(1,1)*T(5,1)+T(5,2)*T(1,2)) /CFAC 

TNOR (6,1)=(T(1,1)*T(6,1)+T(6,2)*T(1,2)) /CFAC 

TNOR(l,2)=0.D0 

TNOR(2,2)=(T(2,2)*T(l,l)-T(2,l)*T(l,2) ) /CFAC 

TNOR(3,2)=(T(3,2) *T(1# 1) -T (3 , 1) *T(1, 2 ) ) /CFAC 

TNOR(4,2)=(T(4,2) *T ( 1 , 1 ) -T (4 , 1 ) *T ( 1 , 2 ) ) /CFAC 

TNOR (5,2)=(T(5,2)*T(1,1)-T(5,1)*T(1,2)) /CFAC 

TNOR(6,2)=(T(6,2)*T(l,l)-T(6,l)*T(l,2) ) /CFAC 


DO  1  1  =  1,  6 

DO  2  J=l,2 

T(I, J)=TNOR(I, J) 

CONTINUE 
CONTINUE 
RETURN 
END 


SUBROUTINE  MULT (TINV, A, T, A2 ) 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
DIMENSION  TINV  (6,6),A(6,6),T(6,6),A1(6,6),A2(6,6) 
DO  1  1=1,6 
DO  2  J=l, 6 

A1(I, J)=0.D0 
A2 (I, J)=0.D0 
CONTINUE 
CONTINUE 
DO  3  1=1, 6 
DO  4  J=l, 6 
DO  5  K=l,6 

A1(I, J)=A(I,K) *T(K, J)+A1(I,  J) 
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5       CONTINUE 
4     CONTINUE 
3      CONTINUE 

DO  6  1=1, 6 
DO  7  J=l, 6 
DO  8  K=l, 6 

A2  (I, J)=TINV(I,K)  *A1  (K,  J) +A2  (I,  J! 
8       CONTINUE 


7     CONTINUE 
6      CONTINUE 

DO  11  1=1,6 
C     WRITE  (*,101)  (A(I,J) ,J=1,6) 

11  CONTINUE 

DO  12  1=1,6 
C     WRITE  (*,101)  (T(I, J) , J=l,6) 

12  CONTINUE 

DO  10  1=1,6 
C     WRITE  (*,101)  (A2 (I, J) , J=l,6) 
1 0     CONTINUE 
C     WRITE  (*,101)  A2(l,l) 

RETURN 
101    FORMAT  (4E15.5) 

END 
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APPENDIX  B 
CRITICAL  VALUE  OF  C  FOR  [X,  X]  BASIS 


PROGRAM  NCCRIT 
HOPF  BIFURCATIONS 
NOMOTO'S  FIRST  ORDER  MODEL 


CALCULATIONS  FOR   CCRITICAL 


C 
C 
C 
C 
C 
C234567 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
DOUBLE  PRECISION  Kl , K2 , K, LPSI, LY, LR, IZ, L, 
&        MASS , NV, NR, NVDOT, NRDOT, NDRS , NDRB , KPSI 


C 


KR,KY,K3 


DIMENSION  AMAT (6,6), FV1 ( 6 ) , IV1 ( 6 ) , YYY (6,6) 
DIMENSION  WR  (  6  )  ,  WI ( 6 )  , TSAVE (6,6), TLUD (6,6), IVLUD ( 6 ) 
DIMENSION  ASAVE ( 6 , 6 ) , Al ( 6 , 6 ) , A2 ( 6 , 6 ) 
OPEN  (11,FILE='CVWN1.RES' ,STATUS='NEW' ) 
FILE= ' CVWN2 . RES ' , STATUS= ' NEW ' ) 
FILE= ' CVWN3 . RES ' , STATUS= ' NEW ' ) 


SVLUD ( 6 ) 


(12 
(13 


OPEN 

OPEN 

WRITE 

READ 

WRITE 

READ 

WRITE 

READ 

WEIGHT=43  5 

IZ 

L 

RHO 

G 

XG 

MASS 

MASS 


ENTER  MIN,MAX,  AND  INCREMENTS  IN  CCRIT' 

CMIN,CMAX,IC 

ENTER  MIN,MAX,  AND  INCREMENTS  IN  WN' 

WNMIN,WNMAX, INCR 

ENTER  WNO' 

WNO 


0 


=  45.0 

=7.3 

=  1.94 

=  32.2 

=0.0104 

=WEIGHT/G 

=MASS/ (0.5*RHO*L**3) 
IZ     =IZ/ (0.5*RHO*L**5) 
XG     =XG/L 
YRDOT  =-0.00000 
YVDOT  =-0.03430 
YR     =+0.00000 
YV     =-0.10700 
YDRS   =+0.01241 
YDRB   =+0.01241 
NRDOT  =-0.00047 
NVDOT  =-0.00000 
NR     =-0.00390 
NV     =-0.00000 
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NDRS   =-0.337*YDRS 
NDRB   =+0.2  83*YDRS 
DH   = (IZ-NRDOT) * (MASS-YVDOT) - 
&      (MASS*XG-YRDOT) * (MASS*XG-NVDOT) 
All= ( (IZ-NRDOT) *YV- (MASS*XG-YRDOT) *NV) /DH 


A12=( (IZ-NRDOT) * (-MASS+YR) - 
&      (MASS*XG-YRDOT) * ( -MASS*XG+NR) ) /DH 

A21  = ( (MASS-YVDOT) *NV- (MASS*XG-NVDOT) *YV) /DH 

A22=( (MASS-YVDOT) * ( -MASS*XG+NR) - 
&      (MASS*XG-NVDOT) * ( -MASS+YR) ) /DH 

Bll= ( (IZ-NRDOT) *YDRS- (MASS*XG-YRDOT) *NDRS) /DH 

B12  = ( (IZ-NRDOT) *YDRB- (MASS*XG-YRDOT) *NDRB) /DH 

B21= ( (MASS-YVDOT) *NDRS- (MASS*XG-NVDOT) *YDRS) /DH 

B22= ( (MASS-YVDOT) *NDRB- (MASS*XG-NVDOT) *YDRB) /DH 

Bl  =B11-B12 

B2  =B21-B22 

Cl= (A11*A22-A21*A12) * (A21*B1-A11*B2 ) 

C2= (A11+A22 ) * (A21*B1-A11*B2 ) +B2* (A11*A22-A21*A12 

C3=- (A21*B1-A11*B2) **2 

A=C1/C2 

B=C3/C2 
C 

EPS=l.D-5 

ILMAX=1500 
C 

DO  1  11=1, INCR 
C 

WN   =WNMIN+ (WNMAX-WNMIN) * (II-l) / (INCR-1) 
C      print  *,wn 

ALPHA0=WN**3 

ALPHA1=2.15*WN**2 

ALPHA2=1.75*WN 

KPSI=-ALPHA1/B 

KY   =-ALPHA0/B 

KR   =- (ALPHA2+A) /B 

GAMA0=WNO**3 

GAMA1=2 . 15*WNO**2 

GAMA2=1.75*WNO 

LY=A+GAMA2 

LPSI=A*LY+GAMA1 

LR=A*LPSI+GAMAO 
C234567 

DO    2    J=1/IC 

CCRIT=CMIN+(J-1)  *  (CMAX-CMIN)  /  (IC-1) 
C  A  MATRIX 

AMAT(1, 1) =0.0 

AMAT(1,2) =1.0 

AMAT(1,3) =0.0 
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AMATd 

,4 

1=0.0 

AMATd 

,5 

1=0.0 

AMATd 

,6 

1=0.0 

AMAT(2 

,1 

1=0.0 

AMAT(2 

,2 

=A 

AMAT(2 

,3 

=  0.0 

AMAT(2 

4 

=B*CCRIT*KPSI 

AMAT(2 

5 

=B*KR 

AMAT(2 

6 

=B*KY 

AMAT(3 

,1 

=  1.0 

AMAT(3 

2 

=  0.0 

AMAT(3 

,3 

=  0.0 

AMAT(3 

,4 

=  0.0 

AMAT(3 

,5 

=  0.0 

AMAT(3 

,6 

=  0.0 

AMAT(4 

1 

=  0.0 

AMAT(4 

,2 

=  0.0 

AMAT(4 

3 

=LPSI 

AMAT(4 

4 

=  0.0 

AMAT(4 

,5 

=  1.0 

AMAT(4 

6 

=-LPSI 

AMAT(5 

1 

=  0.0 

AMAT(5 

2 

=  0.0 

AMAT(5 

3 

=LR 

AMAT(5 

4 

=B*CCRIT*KPSI 

AMAT(5 

5 

=B*KR 

AMAT(5 

6 

=-LR+B*KY 

AMAT(6 

1 

=  0.0 

AMAT(6, 

2 

=  0.0 

AMAT(6, 

3 

=LY 

AMAT(6, 

4' 

=  1.0 

AMAT(6, 

5; 

=  0.0 

AMAT(6, 

6] 

=  -LY 

CALL  RG(6, 6,AMAT,WR/WI/ 0, ZZZ, IV1 , FV1 , IERR) 

CALL  DSTABL(DEOS,WR,WI,FREQ) 

U=CCRIT 

IF  (J.GT.l)  GO  TO  10 

DEOSOO=DEOS 

UOO=U 

LL=0 

GO  TO  2 

DEOSNN=DEOS 

UNN=U 

PR=DEOSNN*DEOSOO 

IF  (PR.GT.0.D0)  GO  TO  3 

LL=LL+1 

IF  (LL.GT.3)  STOP  1000 

IL=0 

UO=UOO 
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UN=UNN 

DEOSO=DEOSOO 

DEOSN=DEOSNN 

UL=UO 

UR=UN 

DEOSL=DEOSO 

DEOSR=DEOSN 

U=(UL+UR) /2.D0 

CCRIT=U 

AMAT (1,1 

(=0.0 

AMATd 

2 

=1.0 

AMAT(1 

3 

=  0.0 

AMAT(1 

4 

1=0.0 

AMATd 

5 

1=0.0 

AMATd 

6 

1=0.0 

AMAT (2 

1 

1=0.0 

AMAT (2 

2 

l=A 

AMAT (2 

3 

1=0.0 

AMAT (2 

4 

l=B*CCRIT*KPSI 

AMAT (2 

5 

1 =B*KR 

AMAT (2 

6 

1 =B*KY 

AMAT (3 

1 

1=1.0 

AMAT (3 

2 

1=0.0 

AMAT (3 

3 

1=0.0 

AMAT (3 

,4 

1=0.0 

AMAT (3 

5 

1=0.0 

AMAT (3 

6 

1=0.0 

AMAT (4 

1 

1=0.0 

AMAT (4 

2 

1=0.0 

AMAT (4 

3 

l=LPSI 

AMAT (4 

4 

1=0.0 

AMAT (4 

5 

1=1.0 

AMAT (4 

6 

>=-LPSI 

AMAT (5, 

1 

=  0.0 

AMAT (5, 

2 

=  0.0 

AMAT (5, 

3 

=LR 

AMAT (5 

4 

=B*CCRIT*KPSI 

AMAT(5, 

5 

=B*KR 

AMAT(5, 

6 

=-LR+B*KY 

AMAT (6, 

1 

=  0.0 

AMAT (6, 

2 

=  0.0 

AMAT(6, 

3 

=LY 

AMAT (6, 

4 

=  1.0 

AMAT(6, 

5 

=  0.0 

AMAT(6, 

6 

=  -LY 

CALL  RG(6, 6,AMAT,WR,WI, 0 ,  ZZZ, IV1 , FV1 , IERR! 

CALL    DSTABL(DEOS,WR,WI,FREQ) 

U=CCRIT 

DEOSM=DEOS 
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UM=U 

PRL=DEOSL*DEOSM 

PRR=DEOSR*DEOSM 

IF  (PRL.GT.O.DO)  GO  TO  5 

UO=UL 

UN=UM 

DEOSO=DEOSL 

DEOSN=DEOSM 

IL=IL+1 

IF  (IL.GT.ILMAX)  STOP  3100 

DIF=DABS(UL-UM) 

IF  (DIF.GT.EPS)  GO  TO  6 

U=UM 

GO  TO  4 

IF  (PRR.GT.0.D0)  STOP  32  00 

UO=UM 

UN=UR 

DEOSO=DEOSM 

DEOSN=DEOSR 

IL=IL+1 

IF  (IL.GT.ILMAX)  STOP  3100 

DIF=DABS(UM-UR) 

IF  (DIF.GT.EPS)  GO  TO  6 

U=UM 

LLL=10+LL 

CCRIT=U 

WRITE  (LLL,*)  CCRIT,WN 

UOO=UNN 

DEOSOO=DEOSNN 

CONTINUE 

CONTINUE 

FORMAT  (215) 
END 


SUBROUTINE  DSTABL ( DEOS , WR , WI , OMEGA ) 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
DIMENSION  WR ( 6 )  ,  WI  ( 6 ) 
DEOS=-1.0D+20 
DO  1  1=1,6 

IF  (WR(I) .LT.DEOS)  GO  TO  1 

DEOS=WR(I) 

IJ=I 
CONTINUE 
OMEGA=WI (IJ) 
OMEGA=DABS (OMEGA) 
RETURN 
END 
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APPENDIX  C 
HOPF  BIFURCATION  PROGRAM  FOR  [X,X]  BASIS 


PROGRAM  KKPSI 
C      HOPF  BIFURCATIONS 
C      NOMOTO'S  FIRST  ORDER  MODEL 
C 

C      CALCULATIONS  FOR  K  AND  CCRITICAL  IF  KPSI  CHANGES  W/  C 
C 
C234567 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
C 

DOUBLE  PRECISION  Kl, K2 , K, LPSI , LY, LR, IZ, L, 
&  MASS , NV, NR, NVDOT, NRDOT, NDRS , NDRB, KPSI ,  KR,  KY,  K3  , 

&  L21  ,L22,L23,L24,L31  ,L32,L33  ,L34,L51  , L52, L53, L54, 

&  Mil ,  M12  ,  M13  ,  M14 ,  M15 ,  M16 ,  M21 ,  M22  ,  M23  ,  M24 ,  M25 ,  M2 6 , 

&  M31 , M32 , M33  , M34 , M35 , M36 , M41 , M42 , M43  , M44 , M45 , M46 , 

&  M51 ,  M52 , M53  ,  M54  , M55 , M56 , M61 ,  M62 ,  M63  ,  M64 , M65 ,  M66 , 

&  Nil , N12 , N13 , N14 , N15 , N16 , N21 , N22 , N23 , N24 , N25 , N2  6 , 

&  N31 , N32 , N33 , N34 , N35 , N3  6 , N41 , N42 , N43 , N44 , N45 , N46 , 

Sc  N51 ,  N52  ,  N53 ,  N54 ,  N55 ,  N5 6 ,  N61 ,  N62  , N63  ,  N64 ,  N6 5 ,  N6 6 

C 

DIMENSION  AMAT  (6,6)^(6,6),  TINV  (6,6),  FV1  (  6 )  ,  IV1  (  6 )  ,  YYY  (6,6) 
DIMENSION  WR  (  6  )  ,  WI  ( 6  )  ,  TSAVE  (6,6),  TLUD  (6,6),  IVLUD  (  6  )  ,  SVLUD  (  6 ) 
DIMENSION  ASAVE ( 6 , 6 ) , Al ( 6 , 6 ) , A2 ( 6 , 6 ) 
OPEN  (10,FILE=/CVWN1.RES' , STATUS= 'OLD' ) 
OPEN  (11,FILE='AKPSI.MAT' ,STATUS='NEW ) 
C       OPEN  ( 1 2, FILE='IREAL. MAT ', STATUS =' NEW ) 
C       OPEN  ( 1 3, FILE='RVALS. MAT' , STATUS =' NEW ) 
C       OPEN  (14, FILE='KKKKK. MAT ', STATUS =' NEW ) 
WEIGHT=435.0 


IZ 

=  45.0 

L 

=7.3 

RHO 

=  1.94 

G 

=32.2 

:g 

=0.0104 

i-iASS 

=WEIGHT/G 

MASS 

=MASS/ (0.5*RHO*L**3) 

IZ 

=IZ/ (0.5*RHO*L**5) 

XG 

=XG/L 

YRDOT 

=-0.00000 

YVDOT 

=-0.03430 

YR 

=+0.00000 

YV 

=-0.10700 

YDRS 

=+0.01241 

YDRB 

=+0.01241 
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NRDOT  =-0.00047 

NVDOT  =-0.00000 

NR     =-0.00390 

NV     =-0.00000 

NDRS   =-0.337*YDRS 

NDRB   =+0.283*YDRS 

DH   = ( IZ-NRDOT) * (MASS-YVDOT) - 
&      (MASS*XG-YRDOT) * ( MAS S*XG -NVDOT) 

All= ( (IZ-NRDOT) *YV- (MASS*XG-YRDOT) *NV) /DH 

A12=( (IZ-NRDOT) * (-MASS+YR) - 
&      (MASS*XG-YRDOT) * ( -MASS*XG+NR) ) /DH 

A21= ( (MASS-YVDOT) *NV- (MASS *XG- NVDOT) *YV) /DH 

A22= ( (MASS-YVDOT) * ( -MASS*XG+NR) - 
&      ( MAS S*XG -NVDOT) * ( -MASS+YR) ) /DH 

Bll=( (IZ-NRDOT) *YDRS-(MASS*XG-YRDOT) *NDRS) /DH 

B12=( (IZ-NRDOT) *YDRB- (MASS*XG-YRDOT) *NDRB) /DH 

B21= ( (MASS-YVDOT) *NDRS- ( MAS S*XG -NVDOT) *YDRS) /DH 

B22= ( (MASS-YVDOT) *NDRB- ( MAS S*XG -NVDOT) *YDRB) /DH 

Bl  =B11-B12 

B2  =B21-B22 
C 

200  WRITE  (*,1004) 

READ  [*,*)       IWN 

INCR=IWN 

WRITE  (*,*)  'ENTER  OBSERVER  WN' 

READ  ( * , * )  WNO 
205  WRITE(*,1007) 

READ  ( * , * )  A3 
C     DO  is  Dsat 

50  WRITE  (*,1006) 

READ   ( * , * )     DO 

Cl= (A11*A22-A21*A12) * (A21*B1-A11*B2 ) 

C2=(A11+A22) *(A21*B1-A11*B2)+B2*(A11*A22-A21*A12 

C3=- (A21*B1-A11*B2) **2 

A=C1/C2 

B=C3/C2 
C      A3 =0.0 
C 
C 


C 
C 


DO  1  11=1, INCR 


READ  (10,*)  CCRIT,WN 

print  *,wn 
ALPHA0=WN**3 
ALPHA1=2 . 15*WN**2 
ALPHA2=1.75*WN 
KPSI=-ALPHA1/B 
KY   =-ALPHA0/B 


48 


KR   =- (ALPHA2+A) /B 
GAMA0=WNO**3 
GAMA1=2 . 15*WNO**2 
GAMA2=1.75*WNO 
LY=A+GAMA2 
LPSI=A*LY+GAMA1 
LR=A*LPSI+GAMA0 
C234567 


k.     MATRIX 

AMAT 

(1 

,1 

)=0 

.0 

AMAT 

(1 

2 

)=1 

.0 

AMAT 

CI 

,3 

1=0 

.0 

AMAT 

(1 

,4 

1=0 

.0 

AMAT 

(1 

,5 

)=0 

.0 

AMAT 

(1 

6 

1=0 

.0 

AMAT 

(2 

-1 

1=0 

.0 

AMAT 

(2 

,2 

)=A 

AMAT 

[2 

3 

1=0 

.0 

AMAT 

(2, 

4 

1 =B*CCRIT*KPSI 

AMAT 

[2, 

5 

l=B*KR 

AMAT 

[2 

1 =B*KY 

AMAT 

^ 

=  1 

.0 

AMAT 

[3 

w 

i=0 

.0 

AMAT 

[3J 

3 

1=0 

.0 

AMAT 

[3, 

4 

1=0 

.0 

AMAT 

13, 

5 

1=0 

.0 

AMAT 

[3, 

6 

1=0 

.0 

AMAT 

14, 

1 

1=0 

.0 

AMAT 

'A, 

2 

1=0 

.0 

AMAT 

r4, 

3 

=LPSI 

AMAT 

4, 

4 

=  0 

0 

AMAT 

4, 

5 

=  1 

0 

AMAT 

4, 

6 

=-LPSI 

AMAT! 

5, 

1 

=  0 

0 

AMATI 

5, 

2 

=  0 

0 

AMATI 

5, 

3] 

=LR 

AMATI 

5, 

4] 

=B*CCRIT*KPSI 

AMATI 

5, 

5! 

=B*KR 

AMATI 

5, 

6] 

=-LR+B*KY 

AMATI 

6, 

11 

=  0. 

0 

AMATI 

6, 

2) 

=  0. 

0 

AMATI 

6, 

3] 

=LY 

AMATI 

6, 

4! 

=  1. 

0 

AMATI 

6, 

5] 

=  0. 

0 

AMATI 

6, 

6) 

=  -LY 

DO  1] 

.  i=: 

L,6 

DO 

12 

!  C 

1=1, 

6 

I 

k.SAVI 

2(1, 

J ) = AMAT ( I , J 

COI> 

ITINt 

JE 

CONT] 

:ni 

JE 
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CALL  RG(6, 6,AMAT,WR/WI,1/YYY,IV1/FV1/IERR: 

CALL  DSOMEG ( IEV, WR, WI , OMEGA, CHECK) 
C       WRITE  (*,*)  IEV 

WRITE  (12,*)   (WR(IREAL) ,IREAL=1,6) 

OMEGAO  =OMEGA 

DO  5  1=1, 6 

T(I,1)=YYY(I,IEV) 
T(I,2)=-YYY(I,IEV+1) 
5      CONTINUE 

IF(IEV.EQ.l)  GO  TO  13 

IF(IEV.EQ.2)  GO  TO  14 

IF(IEV.EQ.3)  GO  TO  15 

IF(IEV.EQ.4)  GO  TO  16 

IF(IEV.EQ.5)  GO  TO  17 

STOP  3004 

13  DO  21  1=1,6 

T ( I , 3 ) =YYY (1,3) 
T ( I , 4 ) =YYY (1,4) 
T ( I , 5 ) =YYY (1,5) 
T ( I ,  6 ) =YYY (1,6) 

2 1  CONTINUE 
GO  TO  30 

14  DO  22  1=1, 6 

T ( I , 3 ) =YYY (1,1) 
T ( I , 4 ) =YYY (1,4) 
T(I,5)=YYY(I,5) 
T ( I , 6 ) =YYY (1,6) 

22  CONTINUE 
GO  TO  3  0 

15  DO  23  1=1,6 


T ( 1 ,  3 ) =YYY ( I , 

1) 

T ( I ,  4 ) =YYY ( I , 

2) 

T ( I , 5 ) =YYY ( I , 

5) 

T ( I , 6 ) =YYY ( I ( 

6) 

23 

CONTINUE 
GO  TO  3  0 

16 

DO  24  1=1, 6 

T ( 1 ,  3 ) =YYY ( I ( 

1) 

T ( I , 4 ) =YYY ( I , 

2) 

T ( I , 5 ) =YYY ( I ( 

3) 

T ( I , 6 ) =YYY ( I , 

■6) 

24 

CONTINUE 
GO  TO  3  0 

17 

DO  25  1=1, 6 

T ( I , 3 ) =YYY ( I , 

1) 

T ( 1 ,  4 ) =YYY ( I , 

2) 

T ( I , 5 ) =YYY ( I , 

3) 

T ( 1 ,  6 ) =YYY ( I , 

4) 

25 

CONTINUE 

30 

CONTINUE 
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NORMALIZATION  OF  THE  CRITICAL  EIGENVECTOR 

CALL  NORMAL (T) 

INVERT  TRANSFORMATION  MATRIX 


DO  2  1=1, 6 

DO  3  J=l, 6 

TINV(I, J) 

=  0.0 

TSAVE(I, J)=T(I, J) 

CONTINUE 

CONTINUE 

CALL  DLUD (6,6 

, TS AVE , 6 

!,TLUD,  IVLUD) 

DO  4  1=1,6 

IF  ( IVLUD ( I ) . EQ . 0 ) 

STOP 

3003 

CONTINUE 

CALL  DILU (6,6 

,TLUD, IVLUD, 

SVLUD) 

DO  8  1=1, 6 

DO  9  J=l,  6 

TINV(I, J) 

=TLUD ( I , 

J) 

CONTINUE 

CONTINUE 

CHECK   Inv(T) 

*A*T 

IMULT=1 

IF  (IMULT.EQ. 

1 )  CALL 

MULT(TINV,A£ 

IF  (IMULT.EQ. 

0 )  STOP 

P=A2 (3,3) 

Q=A2(4,4) 

R=A2 (5,5) 

S=A2 (6,6) 

WRITE (21, *) 

P,Q,R,S 

DEFINITION  OF  Nij 

N11=TINV(1,1) 
N21=TINV(2,1) 
N31=TINV(3,1) 
N41=TINV(4,1) 
N51=TINV(5,1) 
N61=TINV(6,1) 
N12=TINV(1,2) 
N22=TINV(2,2) 
N32=TINV(3,2) 
N42=TINV(4,2) 
N52=TINV(5,2) 
N62=TINV(6,2) 
N13=TINV(1,3) 
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N23=TINV(2, 

3) 

N33=TINV(3, 

3) 

N43=TINV(4/ 

3) 

N53=TINV(5, 

3) 

N63=TINV(6, 

3) 

N14=TINV(1, 

4) 

N24=TINV(2, 

4) 

N34=TINV(3, 

4) 

N44=TINV(4, 

4) 

N54=TINV(5, 

4) 

N64=TINV(6, 

4) 

N15=TINV(1, 

5) 

N25=TINV(2/ 

5) 

N35=TINV(3, 

5) 

N45=TINV(4/ 

5) 

N55=TINV(5, 

5) 

N65=TINV(6, 

5) 

N16=TINV(1, 

6) 

N26=TINV(2, 

6) 

N36=TINV(3, 

6) 

N46=TINV(4, 

6) 

N56=TINV(5, 

6) 

N66=TINV(6, 

6) 

c 

c 

DEFINITION 

OF 

c 

M11=T(1/1) 
M21=T(2,1) 
M31=T(3,1) 
M41=T(4,1) 
M51=T(5,1) 
M61=T(6,1) 
M12=T(1,2) 
M22=T(2,2) 
M32=T(3,2) 
M42=T(4,2) 
M52=T(5,2) 
M62=T(6,2) 
M13=T(l/3) 
M23=T(2,3) 
M33=T(3,3) 
M43=T(4,3) 
M53=T(5,3) 
M63=T(6,3) 
M14=T(1,4) 
M24=T(2,4) 
M34=T(3,4) 
M44=T(4,4) 
M54=T(5,4) 
M64=T(6,4) 

52 


M15  = 

=T(1, 

5) 

M25: 

=T(2, 

5) 

M35: 

=T(3, 

5) 

M45  = 

=T(4, 

5) 

M55  = 

=T(5, 

5) 

M65  = 

=T(6, 

5) 

M16  = 

=T(1, 

6) 

M2  6  = 

=T(2, 

6) 

M36: 

=T(3, 

6) 

M46: 

=T(4, 

6) 

M56  = 

=T  ( 5 , 

6) 

M66: 

=Ti6, 

6) 

Kl= 

=l./8.* 

K2: 

=3.*A3- 

c 

C       Kl=l./8.* ( (ALPHA2**3)+ALPHA0) / (ALPHA2) 

C       K2=3.*A3-.5* (ALPHA2**2) /ALPHAO 

C       K3=l./ ( (B**2) * (D0**2) ) * (ALPHA2+A) * ( ( ALPHAO /ALPHA2 )+ (A* *2 ) ) 

C       K=K1*(K2+K3) 

C 

C      print  *,  wn,k,ccrit 

C 

BL=B/ (3*D0**2) 
C 

C2345678901234567890123456789012345678901234567890123456789012345 
L21=A3*M21**3-BLMCCRIT**3*KPSI**3*M41**3+KR**3*M51**3  + 

&  KY**3*M61**3+ 

&  3*CCRIT**2*KPSI**2*KR*M41**2*M51+ 

&  3*CCRIT**2*KPSI**2*KY*M41**2*M61+ 

&  3*CCRIT*KPSI*KR**2*M51**2*M41+3*CCRIT*KPSI*KY**2*M61**2*M41+ 

Sc  3*KR*KY**2*M61**2*M51  +  3*KR**2*KY*M51**2*M61  + 

&  6*CCRIT*KPSI*KR*KY*M41*M51*M61) 

L22=A3*M22**3-BL* (CCRIT**3*KPSI**3*M42**3+KR**3*M52**3+ 

&  KY**3*M62**3+ 

&  3*CCRIT**2*KPSI**2*KR*M42**2*M52+ 

&  3*CCRIT**2*KPSI**2*KY*M42**2*M62+ 

Sc  3*CCRIT*KPSI*KR**2*M52**2*M42+3*CCRIT*KPSI*KY**2*M62**2*M42+ 

Sc  3*KR*KY**2*M62**2*M52  +  3*KR**2*KY*M52**2*M62  + 

Sc  6*CCRIT*KPSI*KR*KY*M42*M52*M62) 

L23=3*A3*M21**2*M22-BL*(3*CCRIT**3*KPSI**3*M41**2*M42+ 

&  3*KR**3*M51**2*M52+ 

Sc  3*KY**3*M61**2*M62  + 

&  3*CCRIT**2*KPSI**2*KR* (M41**2*M52+2*M41*M42*M51 ) + 

Sc  3*CCRIT**2*KPSI**2*KY*  (M41**2*M62+2*M41*M42*M61 ) + 

Sc  3*CCRIT*KPSI*KR**2*  (M51**2*M42+2*M51*M52*M41 ) + 

&  3*CCRIT*KPSI*KY**2* (M61**2*M42+2*M61*M62*M41 ) + 

Sc  3*KR*KY**2*  (M61**2*M52+2*M61*M62*M51)  + 

Sc  3*KR**2*KY*  (M51**2*M62+2*M51*M52*M61)  + 

Sc         6*CCRIT*KPSI**KR*KY*  (M41*M51*M62+  (M41*M52+M42*M51 )  *M61)  ) 
L24=3*A3*M21*M22**2-BL* (3*CCRIT**3*KPSI**3*M41*M42**2+ 

Sc  3*KR**3*M51*M52**2  + 

&  3*KY**3*M61*M62**2+ 
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&       3*CCRIT**2*KPSI**2*KR* (M42**2*M51+2*M41*M42*M52 ) + 

&     3*CCRIT**2*KPSI**2*KY* (M42**2*M61+2*M41*M42*M62 ) + 

&     3*CCRIT*KPSI*KR**2* (M52**2*M41+2*M51*M52*M42 ) + 

&       3*CCRIT*KPSI*KY**2* (M62**2*M41+2*M61*M62*M42 ) + 

&     3*KR*KY**2* (M62**2*M51+2*M61*M62*M52) + 

&     3*KR**2*KY* (M52**2*M61+2*M51*M52*M62) + 

&     6*CCRIT*KPSI*KR*KY* (M42*M52 *M61+ (M41*M52+M42*M51 ) *M62) ) 

L31=(-l/6) *M11**3 

L32=(-1/6)*M12**3 

L33=(-l/2) *M11**2*M12 

L34=(-1/2)*M11*M12**2 
C2 3 456789 0123456789 0123456789 012345678901234567 89012345678 9 0123 45 

L51=-BL* (CCRIT**3*KPSI**3*M41**3+KR**3*M51**3+KY**3*M61**3+ 
&       3*CCRIT**2*KPSI**2*KR*M41**2*M51+ 
8c  3*CCRIT**2*KPSI**2*KY*M41**2*M61+ 

&  3*CCRIT*KPSI*KR**2*M51**2*M41+3*CCRIT*KPSI*KY**2*M61**2*M41+ 
Sc  3*KR*KY**2*M61**2*M51+3*KR**2*KY*M51**2*M61+ 

Sc  6*CCRIT*KPSI*KR*KY*M41*M51*M61) 

L52=-BL* (CCRIT**3*KPSI**3*M42**3+KR**3*M52**3+KY**3*M62**3+ 
&       3*CCRIT**2*KPSI**2*KR*M42**2*M52+ 
&       3*CCRIT**2*KPSI**2*KY*M42**2*M62+ 

&  3*CCRIT*KPSI*KR**2*M52**2*M42+3*CCRIT*KPSI*KY**2*M62**2*M42+ 
Sc  3*KR*KY**2*M62**2*M52  +  3*KR**2*KY*M52**2*M62  + 

&     6*CCRIT*KPSI*KR*KY*M42*M52*M62 ) 

L53=-BL* (3*CCRIT**3*KPSI**3*M41**2*M42+3*KR**3*M51**2*M52+ 
Sc  3*KY**3*M61**2*M62  + 

&  3*CCRIT**2*KPSI**2*KR* (M41**2*M52+2*M41*M42*M51 ) + 

Sc  3*CCRIT**2*KPSI**2*KY*(M41**2*M62+2*M41*M42*M61)+ 

Sc  3*CCRIT*KPSI*KR**2*  (M51**2*M42+2*M51*M52*M41 )  + 

&  3*CCRIT*KPSI*KY**2* (M61**2*M42+2*M61*M62*M41 ) + 

Sc  3*KR*KY**2*  (M61**2*M52+2*M61*M62*M51)  + 

Sc  3*KR**2*KY*  (M51**2*M62+2*M51*M52*M61)  + 

Sc  6*CCRIT*KPSI*KR*KY*  (M41*M51*M62+  (M41*M52+M42*M51 )  *M61)  ) 

L54=-BL* (3*CCRIT**3*KPSI**3*M41*M42**2+3*KR**3*M51*M52**2+ 
Sc  3*KY**3*M61*M62**2+ 

Sc  3*CCRIT**2*KPSI**2*KR*  (M42**2*M51  +  2*M41*M42*M52  )  + 

Sc  3*CCRIT**2*KPSI**2*KY*  (M42**2*M61  +  2*M41*M42*M62  )  + 

Sc  3*CCRIT*KPSI*KR**2*  (M52**2*M41+2*M51*M52*M42  )  + 

Sc  3*CCRIT*KPSI*KY**2*  (M62**2*M41+2*M61*M62*M42  )  + 

Sc  3*KR*KY**2*  (M62**2*M51+2+M61*M62*M52)  + 

&  3*KR**2*KY* (M52**2*M61+2*M51*M52*M62 ) + 

Sc  6*CCRIT*KPSI*KR*KY*  (M42*M52*M61+  (M41*M52+M42*M51 )  *M62)  ) 

R11=N12*L21+N13*L31+N15*L51 

R12=N12*L23+N13*L33+N15*L53 

R13=N12*L24+N13*L34+N15*L54 

R14=N12*L22+N13*L32+N15*L52 

R21=N22*L21+N23*L31+N2  5*L51 

R22=N22*L2  3+N2  3*L3  3+N25*L53 

R23=N22*L24+N23*L34+N2  5*L54 

R24=N22*L22+N23*L32+N2  5*L52 
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WRITE  (13,*)  Rll  ,R13,R22,R24 
K=(3*R11+R13+R22+3*R24) /8 
WRITE  (11,2001)  WN,K,CCRIT 
1  CONTINUE 
STOP 
1004  FORMAT  ('  ENTER  NUMBER  OF  DATA') 

1006  FORMAT  ('  ENTER  DELTAS AT.') 

1007  FORMAT  ('ENTER  A3') 
2001  FORMAT  (3E15.5) 

END 


SUBROUTINE  DSOMEG ( I JK , WR , WI , OMEGA , CHECK ; 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
DIMENSION  WR ( 6 ) , WI ( 4 ) 
CHECK=-1.0D+25 
DO  1  1=1,6 

IF  (WR(I) .LT. CHECK)  GO  TO  1 

CHECK=WR(I) 

IJ=I 
CONTINUE 

OMEGA=DABS (WI (IJ) ) 
IF  (WI(IJ) .GT.0.D0)  IJK=IJ 
IF  (WI(IJ) .LT.0.D0)  IJK=IJ-1 
RETURN 
END 


SUBROUTINE  NORMAL (T) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  T ( 6 , 6 ) , TNOR (6,6) 

CFAC=T(1,1) **2+T(l,2) **2 

IF  (DABS(CFAC) . LE . (l.D-10) )  STOP  4001 

TNOR(l,l)=l.D0 

TNOR(2,l)=(T(l,l)*T(2,l)+T(2,2) *T(1,2) ) /CFAC 

TNOR (3,1)=(T(1,1)*T(3,1)+T(3,2)*T(1,2)) /CFAC 

TNOR (4,1)=(T(1,1)*T(4,1)+T(4,2)*T(1,2)) /CFAC 

TNOR (5,1)=(T(1,1)*T(5,1)+T(5,2)*T(1,2)) /CFAC 

TNOR (6,1)=(T(1,1)*T(6,1)+T(6,2)*T(1,2)) /CFAC 

TNOR(l,2)=0.D0 

TNOR(2,2)=(T(2,2)*T(l,l) -T (2 , 1 ) *T ( 1 , 2 ) ) /CFAC 

TNOR(3,2)=(T(3,2)*T(l,l)-T(3,l)*T(l,2) ) /CFAC 

TNOR(4,2)=(T(4,2) *T ( 1 , 1 ) -T (4 , 1 ) *T ( 1 , 2 ) ) /CFAC 

TNOR (5,2)=(T(5,2)*T(1,1)-T(5,1)*T(1,2)) /CFAC 

TNOR (6/2)=(T(6,2)*T(l,l)-T(6,l)*T(l,2)) /CFAC 

DO  1  1=1,6 

DO  2  J=l,2 

T(I, J)=TNOR(I, J) 

CONTINUE 
CONTINUE 
RETURN 
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END 
C====================================================== 

SUBROUTINE  MULT (TINV, A, T, A2 ) 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
DIMENSION  TINV (6, 6),A(6,6),T(6,6), Al (6,6), A2 (6,6 
DO  1  1=1,6 
DO  2  J=l, 6 

Al (I, J) =0.D0 
A2 (I, J)=0.D0 

2  CONTINUE 

I  CONTINUE 
DO  3  1=1,6 

DO  4  J=l, 6 
DO  5  K=l, 6 

A1(I, J)=A(I,K) *T(K, J)+A1(I, J) 

5  CONTINUE 
4     CONTINUE 

3  CONTINUE 
DO  6  1=1,6 

DO  7  J=l,6 
DO  8  K=l, 6 

A2 (I, J)=TINV(I,K) *A1(K, J)+A2 (I, J) 
8       CONTINUE 
7     CONTINUE 

6  CONTINUE 

DO  11  1=1, 6 
C     WRITE  (*,101)  (A(I,J) ,J=1,6) 

II  CONTINUE 

DO  12  1=1,6 
C     WRITE  (*,101)  (T(I, J) ,J=1,6) 
12     CONTINUE 

DO  10  1=1,6 
C     WRITE  (*,101)  (A2 (I, J) ,J=1,6) 
1 0     CONTINUE 
C     WRITE  (*,101)  A2(l,l) 

RETURN 
101    FORMAT  (4E15.5) 

END 
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